#include "simulation.h"
#include <string.h>

void sim_delete( Simulation *sim )
{
    CHECK_SIM(sim);
	field_delete(sim->eta1);
	field_delete(sim->zeta1);
	field_delete(sim->dvg_x);
	field_delete(sim->dvg_y);
	field_delete(sim->dvg_x_bis);
	field_delete(sim->dvg_y_bis);
	field_delete(sim->sigma11); 
	field_delete(sim->sigma12);
	field_delete(sim->sigma22);
	field_delete(sim->sigma11_bis);
	field_delete(sim->sigma12_bis);
	field_delete(sim->sigma22_bis);
	field_delete(sim->ry);
	field_delete(sim->rx);
    field_delete(sim->Lrhu);
	field_delete(sim->Lrhv);
	field_delete(sim->Lrho); 
	field_delete(sim->dvg); 
	field_delete(sim->Jv_pred);
	field_delete(sim->Ju_pred);
    field_delete(sim->Jv);
    field_delete(sim->Ju);
	field_delete(sim->delta_Jv1);
	field_delete(sim->delta_Ju1);
	field_delete(sim->delta_Jv);
	field_delete(sim->delta_Ju);
    field_delete(sim->v);
    field_delete(sim->u);
    field_delete(sim->gradPy);
    field_delete(sim->gradPx);
    field_delete(sim->P);
    field_delete(sim->rho);
    
    memset( sim, 0, sizeof(Simulation) );
    free(sim);
}

Simulation *sim_create(const int    nx,
                       const int    ny,
                       const double Lx,
                       const double Ly,
					   const double dt,
					   const double K, 
					   const double eta,
					   const int ndet_i, 
					   const int ndet_j)
{
    assert(nx>1);
    assert(ny>1);
    assert(Lx>0);
    assert(Ly>0);
    
    {
        //----------------------------------------------------------------------
        // create the simulation struct
        //----------------------------------------------------------------------
        Simulation *sim =
        (Simulation *) calloc( 1, sizeof(Simulation));
        if( NULL == sim )
        {
            fprintf( stderr, "can't allocation Simulation");
            exit(1);
        }
        
        //----------------------------------------------------------------------
        // fill the sim
        //----------------------------------------------------------------------
        sim->nx      = nx;
        sim->ny      = ny;
        sim->Lx      = Lx;
        sim->Ly      = Ly;
        sim->dx      = Lx/((double)(nx));
        sim->dy      = Ly/((double)(ny));
        sim->dx2     = sim->dx * sim->dx;
        sim->dy2     = sim->dy * sim->dy;
        sim->inv_dx  = 1.0 / sim->dx;
        sim->inv_dy  = 1.0 / sim->dy;
        sim->inv_dx2 = 1.0 / sim->dx2;
        sim->inv_dy2 = 1.0 / sim->dy2;
		sim->dt      = dt;
		sim->K       = K; 
		sim->eta     = eta;
		sim->ndet_i  = ndet_i; 
		sim->ndet_j  = ndet_j; 
		
        sim->eq = vdw_fluid;
		sim->visco = eta_1; 
		sim->second_visco = zeta_1; 
        sim->T  = 0.95;
        
        sim->wx = (2*M_PI)/sim->Lx;
        sim->wy = (2*M_PI)/sim->Ly;
        
        //----------------------------------------------------------------------
        // allocate the fields
        //----------------------------------------------------------------------
#define __CREATE(NAME) sim->NAME = field_create( #NAME, sim->nx, sim->ny )
        __CREATE(rho);
        __CREATE(P);
        __CREATE(gradPx);
        __CREATE(gradPy);
        __CREATE(u);
        __CREATE(v);
		__CREATE(delta_Ju);
		__CREATE(delta_Jv);
		__CREATE(delta_Ju1);
		__CREATE(delta_Jv1);
        __CREATE(Ju);
        __CREATE(Jv);
		__CREATE(Ju_pred);
		__CREATE(Jv_pred);
		__CREATE(dvg); 
		__CREATE(Lrho);
		__CREATE(Lrhv); 
        __CREATE(Lrhu);
        __CREATE(rx);
		__CREATE(ry);
		__CREATE(sigma22); 
		__CREATE(sigma12);
		__CREATE(sigma11);
		__CREATE(sigma11_bis);
		__CREATE(sigma12_bis);
		__CREATE(sigma22_bis);
		__CREATE(dvg_y_bis);
		__CREATE(dvg_x_bis);
		__CREATE(dvg_y);
		__CREATE(dvg_x);
		__CREATE(eta1);
		__CREATE(zeta1);
        return sim;
    }
}
